Loading Library

library(raster)
library(horizon)
library(rgdal)
pro=CRS("+init=epsg:28992")
WGS84<-CRS("+init=epsg:4326")

Loading the data

The altitude data has a resolution of 90 meters.

#r <- getData('alt', country='Austria') #used for downloading the data
r <- raster("data/AUT_msk_alt.grd")
plot(r,main="Altitude map of Austria")

print(r)
## class       : RasterLayer 
## dimensions  : 336, 948, 318528  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 9.4, 17.3, 46.3, 49.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : /nobackup/users/dirksen/SVF_highres/SVF/data/AUT_msk_alt.grd 
## names       : AUT_msk_alt 
## values      : 108, 3533  (min, max)

Loading the data

#r.NED <- getData('alt', country='NL')# used for downloading the data
r.NED <- raster("data/NLD_msk_alt.grd")
plot(r.NED,main="Altitude map of Netherlands")

print(r.NED)
## class       : RasterLayer 
## dimensions  : 348, 480, 167040  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 3.3, 7.3, 50.7, 53.6  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : /nobackup/users/dirksen/SVF_highres/SVF/data/NLD_msk_alt.grd 
## names       : NLD_msk_alt 
## values      : -8, 286  (min, max)

Calculating the SVF

s <- svf(r, nAngles=16, maxDist=500, ll=TRUE)
## 0%...10%...15%...20%...25%...30%...35%...40%...45%...50%...55%...60%...65%...70%...75%...80%...85%...90%...95%...
s.NED <- svf(r.NED, nAngles=16, maxDist=500, ll=TRUE)
## 0%...10%...15%...20%...25%...30%...35%...40%...45%...50%...55%...60%...65%...70%...75%...80%...85%...90%...95%...
print(s)
## class       : RasterLayer 
## dimensions  : 336, 948, 318528  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 9.4, 17.3, 46.3, 49.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : /tmp/Rtmp0NECCC/raster/r_tmp_2016-11-24_151201_15395_28883.grd 
## names       : layer 
## values      : 0.7161961, 1  (min, max)
print(s.NED)
## class       : RasterLayer 
## dimensions  : 348, 480, 167040  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 3.3, 7.3, 50.7, 53.6  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : /tmp/Rtmp0NECCC/raster/r_tmp_2016-11-24_151207_15395_49376.grd 
## names       : layer 
## values      : 0.9967596, 1  (min, max)

Plotting routine

plot(s,main="SVF for Austria")

plot(s.NED,main="SVF for Netherlands")

A higher resolution example: AHN

The raw data from the AHN is available with this link. As an example I included the map number r50fn2.

Besides this product from PDOK and nationaal Georegister also a 3D map of the Netherlands is available

xml.link.raw<-"http://geodata.nationaalgeoregister.nl/ahn2/atom/ahn2_05m_ruw.xml" 
r.NED<-raster("AHN/r09dn1.tif") #runnning the 0.5m tif takes approximately 2 hours
plot(r.NED,main="Altitude map of Netherlands")

print(r.NED)
## class       : RasterLayer 
## dimensions  : 12500, 10000, 1.25e+08  (nrow, ncol, ncell)
## resolution  : 0.5, 0.5  (x, y)
## extent      : 110000, 115000, 556250, 562500  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs 
## data source : /nobackup/users/dirksen/SVF_highres/SVF/AHN/r09dn1.tif 
## names       : r09dn1

Regridding the data

The runtime for a 0.5m resolution grid is long. Therefore I regrid the raster to a 5m grid. Calculating the SVF will now take approximately 1 minute.

r.NED<-aggregate(r.NED,fact=10) #regridding to a 5m grid, this takes just a minute to run
res(r.NED)
## [1] 5 5
plot(r.NED)

system.time(
s.NED <- svf(r.NED, nAngles=16, maxDist=100, ll=F)
)
## 0%...10%...15%...20%...25%...30%...35%...40%...45%...50%...55%...60%...65%...70%...75%...80%...85%...90%...95%...
##    user  system elapsed 
## 125.513   0.109 126.084
plot(s.NED,main="SVF for Netherlands")

Other area’s

see for the ahn_units.

ahn.units<-readOGR(dsn="AHN_kaartbladen",layer="ahn_units")
## OGR data source with driver: ESRI Shapefile 
## Source: "AHN_kaartbladen", layer: "ahn_units"
## with 1441 features
## It has 3 fields
proj4string(ahn.units)<-pro
plot(ahn.units)

coords.knmi<-readRDS("data/coordsKNMI.rda")
coords.knmi<-na.omit(coords.knmi)
coords.knmi<-data.frame(coords.knmi)
coordinates(coords.knmi)<-~DS_LON+DS_LAT
proj4string(coords.knmi)<-WGS84
coords.knmi<-spTransform(coords.knmi,pro)


ov<-over(coords.knmi,ahn.units)
combined<-cbind(ov,as.data.frame(coords.knmi))
print(combined[which(combined$DS_CODE=="260_H"),])
##      LO_X   LO_Y  UNIT DS_CODE   DS_LON   DS_LAT
## 10 140000 456250 32cn1   260_H 140783.3 456758.2

For station the bilt

r.Bilt<-raster("AHN/r32cn1.tif") #runnning the 0.5m tif takes approximately 2 hours
plot(r.Bilt,main="Altitude map of De Bilt")
points(coords.knmi,col="red",lwd=3)

print(r.Bilt)
## class       : RasterLayer 
## dimensions  : 12500, 10000, 1.25e+08  (nrow, ncol, ncell)
## resolution  : 0.5, 0.5  (x, y)
## extent      : 140000, 145000, 456250, 462500  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs 
## data source : /nobackup/users/dirksen/SVF_highres/SVF/AHN/r32cn1.tif 
## names       : r32cn1
r.Bilt<-aggregate(r.Bilt,fact=10) #regridding to a 5m grid, this takes just a minute to run
res(r.Bilt)
## [1] 5 5
plot(r.Bilt)
points(coords.knmi,col="red",lwd=3)

system.time(
s.Bilt <- svf(r.Bilt, nAngles=16, maxDist=100, ll=F)
)
## 0%...10%...15%...20%...25%...30%...35%...40%...45%...50%...55%...60%...65%...70%...75%...80%...85%...90%...95%...
##    user  system elapsed 
## 126.430   0.147 127.049
plot(s.Bilt,main="SVF for De Bilt")
points(coords.knmi,col="red",lwd=1)